A Novel Mesh Generation Algorithm for Field-Level Coupled Flow and Geomechanics Simulations

Generating a suitable hexahedral mesh for field-level coupled flow and geomechanics computation is still a challenging task. The lack of an analytical representation of the reservoir geometry makes generating a valid mesh for geomechanics a tedious and time-consuming endeavor. Indeed, usually the reservoir description comes as a static-model, i.e., corner-point geometry mesh, which is appropriate for flow simulations but not necessarily for mechanics. Having different meshes for flow and mechanics can alleviate this constraint but there is an intrinsic interpolation error in the procedure. In order to minimize this error, we generate in the pay-zone region a mechanics mesh by smoothing the flow mesh. This procedure ensures that a valid mesh for finite element purposes is obtained and the projector, which maps pressures from the flow mesh into the generated reference mechanics mesh, is the identity matrix. This resulting pay-zone mesh is propagated into its surroundings by means of elliptical smoothing using linear elasticity. We also show that the same hexahedral mechanic pay-zone’s mesh can be converted to a valid tetrahedral mesh by exploiting every hexahedron in a certain number of tetrahedrons. This latest procedure is attractive to bringing tetrahedral meshes into the picture, which allows treating constraints for the meshing process, i.e., faults and pitchouts, accordingly. Finally, field-scale reservoir compaction and subsidence computations are carried out by using continuous Galerkin finite elements for mechanics, coupled with a slightly compressible single-phase flow simulator in order to demonstrate the applicability of the proposed algorithm. exactly honored. For general reservoir meshes with faults and pitchouts, that claim cannot be made because some hexahedral elements degenerate (due to pitchouts, for instance), and faults must be smoothed out to avoid dealing with contact problems on those discontinuous surfaces. Angus et al. [8] followed a similar approach pointing out that the resulting mechanics mesh must be watertight, meaning that there must be no gaps in the mesh, and each node must be connected to neighboring nodes properly. For them, faults and other subsurface discontinuities pose a major challenge, and thus they must be properly handled when exported from the geological model. On the other hand, Lie et al. [9] recognized the necessity of employing unstructured meshes when dealing with complex geological features and also to overcome limitations in the numerical schemes such as K-orthogonal grids. Mustapha [10] developed a mesh generation tool aimed to tackle complex fractured geological media by using triangular/tetrahedral grids. Also, Juntunen and Wheeler [11] employed modern mesh generation techniques to improve the quality of the resulting mesh while reducing the number of elements and capture the geometry accurately by using orthogonally optimized hexahedral meshes. These approaches based on unstructured and more general hexahedral meshes suggest different meshes for flow and mechanics in the pay zone. The remainder of the paper is organized as follows. Section 2 presents the mathematical model and the governing partial differential equations. Section 3 focuses on the mesh generation process based on the elasticity operator. Section 4 discusses how to attract the mesh towards given features by using a pressure-drop as driving force. In Section 5, we present concrete numerical examples of mesh generation and coupled flow, and geomechanics simulations are conducted to demonstrate the applicability of the methodology proposed in this paper. Finally, Sections 6 and 7 state concluding remarks and future work respectively. 2. MATHEMATICAL MODEL This section discusses the governing equations for linear homogeneous isotropic poroelasticity and their finite element formulation. We omit details for the sake of brevity, a more detailed treatment can be found [4]. We consider a bounded domain 3 , 2 R ⊂ Ω and its boundary Ω ∂ = Γ . Let h ξ be a non-degenerate, quasi-uniform conforming partition of Ω composed of quadrilaterals or hexahedrons. For the elasticity part, we start from the equilibrium equation for a quasi-steady process: (1) where σ is the stress tensor, n̂ is the outer normal unitary vector. The boundary conditions for mechanics can be assumed to be of Dirichlet type on ΓD u , and Neumann type on u N Γ , where the external tractions are prescribed. Hooke's law and Biot's poroelastic theory define σ by [4]: (2) Where C is the elastic moduli of isotropic elasticity, δ is the Kronecker delta, and λ ,μ are the Lamé constants, and Π is the fourth-order identity tensor. The strain tensor ε is defined by: (3) The Lamé constants can be expressed regarding familiar quantities such as Young's Modulus, E , and Poisson ratio, ν [4]: (4) where G is the Shear Modulus [4]. For the flow part, we are interested in the steady-state continuity equation [4], this is: (5) where K is the absolute permeability tensor, μ is the dynamic viscosity, and p is the fluid. The typical boundary conditions for pressure involve Neumann or no-flow namely: We derive weak forms for Eq. (1) and (2) by multiplying by a test function, v ∈H0 1 Ω ( ) , and integrating over the domain and applying the Gauss divergence theorem. We omit details here for the sake of brevity; a more detailed treatment can be found in [1-4]. This leads to our finite element model for steady state linear isotropic poroelasticity: